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We study the effect that the heterogeneity present among the elements of an ensemble 
of coupled excitable neurons have on the collective response of the system to an external 
signal. We have considered two different interaction scenarios, one in which the neurons 
are diffusively coupled and another in which the neurons interact via pulse-like signals. We 
found that the type of interaction between the neurons has a crucial role in determining the 
response of the system to the external modulation. We develop a mean-field theory based 
on an order parameter expansion that quantitatively reproduces the numerical results in the 
case of diffusive coupling. 

I. INTRODUCTION 

Synchronized behavior arising among the constituents of an ensemble is common in nature. 
Examples include the synchronized flashing of fireflies, the blossoming of flowers, cardiac cells 
giving rise to the pacemaker role of the sinoatrial node of the heart and the electrical pulses of 
neurons. This global behavior can originate from a common response to an external stimulus 
or might appear in autonomous, non-forced, systems. The theoretical basis for the understand- 



ing of synchronization in non-forced systems was put forward by Winfree ( H17H ) who showed 
that the interaction -i.e. coupling- between the constituents is an essential ingredient for the ex- 
istence of a synchronized output. The seminal work on coupled oscillators by Kuramoto (01]) 
offered a model case whose solution confirms the basic hypothesis of Winfree: while interaction 
helps towards the achievement of a common behavior, a perfect order can be achieved only in the 
absence of diversity -heterogeneity- among the components of an ensemble. In the Kuramoto 
model, diversity manifests when the oscillators have different natural frequencies, those they dis- 
play when uncoupled from each other. While this result holds for systems that can be described 
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by coupled oscillators, recent results indicate that in other cases diversity among the constituents 
might actually have a positive role in the setting of a resonant behavior with an external sig- 
nal. This was first demonstrated i n (llw\ ) and has been since extended to many other systems 
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19D). In the case of non-forced excitable systems, a unifying 



treatment of the role of noise and diversity has been developed in fl 1 511 - 

Many biological systems, including neurons, display excitability as a response to an external 
perturbation. Excitability is characterized by the existence of a threshold, a largely independent 
response to a suprathreshold input and a refractory time. It is well established that the dynamical 
features of neurons can be described by excitable models: when a neuron is perturbed by a single 
impulse, the neuron can generate a single spike, and when perturbed by a continuum signal, a 
train of spikes with a characteristic firing frequency can appear instead. Although the creation and 
propagation of electrical signals has been thoroughly studied by physiologists for at least a century, 
the most important landmark in this field is due to Hodgkin and Huxley ([jg]), who developed the 
first quantitative model to describe the evolution of the membrane potential in the squid giant axon. 
Because of the central importance of cellular electrical activity in biological systems and because 
this model forms the basis for the study of excitable systems, it remains to this date an important 
model for analysis. Subsequently a simplified version of this model, known as the FitzHugh- 
Nagumo (FHN) model ([3]), that captures many of the qualitative features of the Hodgkin-Huxley 
model was developed. The FHN model has two variables, one fast and one slow. The fast variable 
represents the evolution of the membrane potential and it is known as the excitation variable. The 
slow variable accounts for the K + ionic current and is known as the recovery variable. One virtue 
of this model is that it can be studied using phase-plane methods, because it is a two variable 
system. For instance, the fast variable has a cubic nullcline while the slow variable has a linear 
nullcline and the study of both nullclines and their intersections allow to determine the different 
dynamical regimes of the model. The FHN model can then be extended by coupling the individual 
elements, and provides an interesting model, for example, of coupled neuron populations. 

In the initial studies of the dynamics of the FHN and similar models, the individual units (e.g. 
neurons) were treated as identical. However, it is evident that real populations of neurons display 
a large degree of variability, both in morphology and dynamical activity; that is, there is a diversity 
in the population of these biological units. Then, it is natural to ask what role diversity plays in the 
global dynamical behavior of these systems and a lot of activity has been developed along these 
lines in the recent years. In general, it has been found, as noted above, that diversity can in fact 
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be an important parameter in controlling the dynamics. In particular, it has been shown that both 
excitable and bistable systems can improve their response to an external stimulus if there is an 
adequate degree of diversity in the constituent units (1 14ID. 

In this paper we continue our study of the effect of diversity, investigating in some depth its 
role in the FHN model of excitable systems. We consider in detail a system of many neurons 
coupled either chemically or electrically. We show that in both cases a right amount of diversity 
can indeed enhance the response to a periodic external stimulus and we discuss in detail the dif- 
ference between the two types of coupling. The outline of the paper is as follows: In section|II]we 
present the FHN model and the coupling schemes considered. Next, in section [III] we present the 
main results we have obtained from numerical simulations. Afterwards, in section [IV] we develop 
an approximate theoretical treatment based on an order parameter expansion which allows us to 
obtain a quantitative description of the behavior of the model and compare its predictions with the 
numerical results. Finally, in section|V]we summarize the conclusions. 

II. FITZHUGH-NAGUMO MODEL 

A. Dynamical equations 

Let us consider a system of excitable neurons described by the FHN model. The dynamical 
equations describing the activity of a single neuron are: 

dx 

e— = x(l - x)(x - b) - y + d, (1) 
dt 

d y ^ ™ 

— = x - cy + a, (2) 
dt 

where x(t) and y(i) represent, respectively, the fast membrane potential and the slow potassium 
gating variable of a neuron. We assume that the time scales of these variables are well separated 
by the small parameter e = 0.01. Other parameters are fixed to b = 0.5, c = 4.6 and d = 0.1 (|40), 
while the value of a will fluctuate from one neuron to the other, so reflecting the intrinsic diversity 
in the neuronal ensemble. 

Let us concentrate first on the dynamics of a single neuron as described by Eqs. ([]])-©■ This 
dynamics has a strong dependence on the parameter a. Three different operating regimes are 
identify: for a < -0.09 the system has a stable focus in the right branch of the cubic nullcline 
leading the system to an excitable regime; for -0.09 < a < 0.01 a limit cycle around an unstable 
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focus appears (oscillatory regime) and for 0.01 < a a stable focus appears again, now at the left 
side of the cubic nullcline (excitable regime). Figure[T](a) shows the nullclines f(x) = x{\ - x)(x - 
b) + d and g(x, a) = (x + a)/c of the system in the three operating regimes described above, for 
a = -0.1, 0.0 and 0.06, respectively. In the excitable regime, spikes (also known as pulses), can 
appear as a result of an external perturbation of large enough amplitude. A convenient definition 
is that a spike appears when the membrane potential of the neuron exceeds a certain threshold 
value, e.g. x > 0.5. In the oscillatory regime, spikes appear spontaneously with an intrinsic firing 
frequency v which, as shown in Figure [T](b), does not depend much on the value of a. 
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FIG. 1 : (a) Nullclines of the FHN system for three different values of the parameter a. f(x) (solid line) and 
g(x, a) for a = -0.1 (dotted line); a = 0.0 (dash-dotted line) and a = 0.06 (dashed line), (b) Dependence of 
the firing frequency v on a. 

To illustrate the dynamical behavior of x(i) and y(t), we show in Figure [2] the phase-portrait for 
three different values of a = -0.1, 0.0 and 0.06 corresponding to the three nullclines represented 
in Figured] (a). 



B. Coupling scenarios: electrical versus chemical interaction 

We consider now that we have an ensemble of iV coupled neurons. Each one is described 

by dynamical variables jc,-(0, y;(0» z = 1 iV, obeying the FitzHugh-Nagumo equations. The 

neurons are not isolated from each other, but interact mutually. The full set of equations is now: 

dx- 

e-± = Xi (l - Xi )(Xi - b) - y t + d + I? n (i), (3) 
dt 

dyi 

— = Xi - cyj + a t , (4) 
dt 
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FIG. 2: Phase-space portrait of the FHN system for different values of a. Gray line represents the evolution 
of {x, y}. f(x) (black line) and g(x, a) for (a) a = -0. 1 dotted line (excitable regime); (b) a = 0.0 dash-dotted 
line (oscillatory regime) and (c) a = 0.06 dash line (excitable regime). 

where lf n {t) is the coupling term accounting for all the interactions of neuron i with other neu- 
rons. To take into account the natural diversity of the units, we assume that the parameter a,, that 
controls the degree of excitability of the neuron, varies from neuron to neuron. In particular we 
assign to a, random values following a Gaussian distribution with mean (a,) = a and correlations 
^(a, - a)(aj - a)^ = 8ij(T 2 . We use cr as a measure of the heterogeneity of the system and, in the 
following, we use the value of cr as an indicator of the degree of diversity. If cr = 0, all neurons 
have exactly the same set of parameters, while large values of cr indicate a large dispersion in the 
dynamical properties of individual neurons. 

The most common way of communication between neurons is via chemical synapses, where 
the transmission is carried by an agent called neurotransmitter. In these synapses, neurons are 
separated by a synaptic cleft and the neurotransmitter has to diffuse to reach the post-synaptic 
receptors. In the chemical coupling case the synaptic current is modeled as: 




In this configuration, we consider that neuron i is connected to N c neurons randomly chosen from 
the set of N - I available neurons. Once a connection is established between neuron i and j, we 
assume that the reciprocal connection is also created. Then, the connection fraction of each neuron 
is defined as / = N C /(N - 1). In Eq. (0) K determines the coupling strength and g, ; represents the 
maximum conductance in the synapse between the neurons i and j. For simplicity, we limit our 
study to the homogeneous coupling configuration, where g, ; = 1 if neurons i and j are connected 
and gij = otherwise. The character of the synapse is determined by the synaptic reversal potential 
of the receptor neuron, E s .. An excitatory (resp. inhibitory) synapse is characterized by a value 
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of E* greater (resp. smaller) than the membrane resting potential. We consider E* = 0.7 for the 



neurons (those that project excitatory synapses) in the system as f e = N e /N being N e the number 
of excitatory neurons. 

Finally, rj(t) is a time dependent function representing the fraction of bound receptors and it is 
given by: 



where a = 2.5 and /3 = 3.5 are the rise and decay time constants, respectively. Here t on = 
0.1 represents the time the synaptic connection remains active and t is the time from the spike 
generation in the presynaptic neuron. 

There is another type of synapse where the membranes of the neurons are in close contact 
and thus the transmission of the signal is achieved directly (electrical synapses). In this case of 
electrically-mediated interactions, also known as diffusive coupling, the total synaptic current is 
proportional to the sum of the membrane potential difference between a neuron and its neighbors, 
and it is given by: 



The last ingredient of our model is the presence of an external forcing that acts upon all neurons 
simultaneously. Although our results are very general, for the sake of concreteness we use a 
periodic forcing of amplitude A and period T. More precisely, the dynamical equations under the 
presence of this forcing are modified as: 



excitatory synapses and E* = -2.0 for the inhibitory ones. We also define the fraction of excitatory 




(6) 




(7) 



e— l - = Xi(l - Xi)(Xi -b)-yi + d + I. 



(t), 



(8) 




(9) 



which is the basis of our subsequent analysis. 



III. RESULTS 



We are interested in analysing the response of the global system to the external forcing. We 
will show that its effect can be enhanced under the presence of the right amount of diversity in the 
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set of parameters a,, i.e. a conveniently defined response will reach its maximum amplitude at an 
intermediate value of the root mean square value cr. 

In order to quantify the global response of the system with respect to the diversity, we use the 
spectral amplification factor defined as 



where X(t) = 4 2Zf=i X M is the global, average collective variable of the system and <•) denotes a 
time average. We will analyze separately the cases of electrical and chemical coupling. 



In this subsection we concentrate on the situation in which the neurons are electrically (diffu- 
sively) coupled in a random network, where a neuron i is connected randomly with N c = f(N - 1) 
other neurons. The mean value of the Gaussian distribution of the parameters a, is fixed to a = 0.06 
and the coupling strength to K = 0.6. Figure [3] shows the spectral amplification factor 77 as a func- 
tion of the diversity a for fixed values of the amplitude A = 0.05 and two different values of the 
period T of the external forcing, for an increasing connection fraction /. It can be seen from Fig- 
ure |3 that intermediate values of cr give a maximum response in the spectral amplification factor. 
Moreover, the maximum value shifts slightly to smaller values of cr as the fraction of connected 
neurons / increases. We have also observed that a period T of the external forcing close to the 
inverse of the intrinsic firing frequency of the neurons (v « 0.9, according to Figure Q})) yields the 
largest response. 

In order to further illustrate the response of the system to the external sinusoidal modulation, Fig- 
ure U] shows the raster plot of the ensemble (lower panels) and the time traces of ten randomly 
chosen neurons (upper panels) for different values of the diversity parameter. It can be seen in 
both the top and bottom panels of this figure that an intermediate level of diversity gives a more 
regular behavior than either smaller or larger values of cr. This fact is more evident in the time 
traces where the amplitude of the oscillation varies randomly for large cr values. 




(10) 



1. Electrical coupling 



2. Chemical coupling 



We consider in this subsection the situation in which the units are chemically coupled. Figure 
|5] shows the spectral amplification factor 77 as a function of the diversity cr for fixed values of the 
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FIG. 3: Spectral amplification factor rj as a function of a for an increasing fraction of connected neurons 
/ for two different periods of the external modulation, (a) T - 1.6 and (b) T - 1.11. Other parameters: 
a = 0.06, K = 0.6 and A = 0.05. 
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FIG. 4: Time traces of ten randomly chosen neurons and raster plot (every time a neuron spikes a dot is 
drawn) of the fully connected, / = 1, ensemble in the case of electrical coupling for three different values 
of the diversity parameter: (a) <x = 0.1, (b) <x = 0.4 and (c) cr = 0.9. Other parameters: a = 0.06, K = 0.6, 
A - 0.05 and T = 1.6. 



amplitude A = 0.05 and two different periods of the external modulation when the fraction of 
randomly connected neurons / increases. The coupling strength is fixed to K = 1.5. The fraction 
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FIG. 5: Spectral amplification factor rj as a function of cr for an increasing fraction of connected neurons /. 
The fraction of excitatory neurons was fixed to f e = 0.8. Two different periods of the external modulation 
have been considered: (a) T = 1.6 and (b) T — 1.11. 

of excitatory neurons in the system is set to f e = 0.8. It can be seen from Figure[5]that the spectral 
amplification factor rj increases as / increases, reaching the maximum response at / =* 0.05. 
Interestingly, beyond this value rj does not change significantly, indicating that the response of the 
system does not improve when the percentage of connected neurons is larger than 5%. Or, put 
in another way, with a 5% connectivity, the system already behaves as being fully connected as 
far as the response to the external forcing is concerned. It is also worth noting that the maximum 
response is always at the same value of cr, independent of /. The effect of changing the ratio 
of excitatory/inhibitory synapses in our system is shown in Figure [6] in the globally coupled case 
/ = 1. The spectral amplification factor increases as the fraction of excitatory connections f e 
increases, while the position of the maximum shifts slightly to larger values of cr 0.05. 

Comparing the results from both electrical and chemical coupling schemes, it can be seen that 
the electrical coupling gives a larger value of rj but requires, at the same time, a larger diversity. 
The electrical coupling also exhibits a larger range of diversity values for which the system has 
an optimal response compared with the chemical coupling. In contrast, the optimal response in 
the chemical coupling scheme occurs for small values of the diversity and does not significantly 
change in amplitude and width when the percentage of connected neurons / is increased above 
5%. 



10 




FIG. 6: Spectral amplification factor 77 as a function of o~ for an increasing fraction of excitatory synapses 
f e in the case that neurons are globally coupled, / = 1. Other parameters: a - 0.06, K = 1.5, A = 0.05 and 
T = 1.6. 

IV. ORDER PARAMETER EXPANSION 



It is possible to perform an approximate analysis of the effect of the diversity in the case of dif- 
fusive (electrical) coupling. The analysis allows to gain insight into the amplification mechanism 
by showing how the effective nullclines of the global variable X(t) are modified when varying <x. 
The theoretical analysis is based upon a modification of the so-called order parameter expansion 
developed by Monte and D'Ovidio 0llll.ll2n) along the lines of [8Q. The approximation begins by 
expanding the dynamical variables around their average values X(t) = jj £ ; x; and Y{t) = jj 2,y,- 
as Xj(t) = X(t) + 6 X j(t), yj(t) = Y(t) + 6 y .{t) and the diversity parameter around its mean value 
a.j = a + 6". The validity of this expansion relies on the existence of a coherent behavior by which 
the individual units x; deviate in a small amount 6 X - from the global behavior characterized by the 
global average variable X(t). It also assumes that the deviations 6" are small. We expand equations 
© for ^ and © for ^ up to second order in 6*, 6 y t and 6f; the resulting equations are: 



dx 1 

e-£ = f (X, Y) + f x (X, Y) 8* + f y (X, Y) 8] + -f xx (X, Y) (5*f , 

^ = g (X, Y, a) + g x (X, Y, a) 8J + g a (X, Y, a) 



(ID 
(12) 
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where 

f(x,y) = x(l - x)(x -b)-y + d-Kx, 
g(x,y, a) = x - cy + a + A sin( yn , 
and we used the notation f x to indicate the derivative of / with respect to x and so forth. Note that 
Eq. (fT2l) is exact since it is linear in all the variables. If we average Eq. (fTTI)-(fT2l using <•> = 4 ■ 
we obtain: 

e^- = f{X,Y) + \f xx {X,Y)Q.\ (14) 
at 2 

^ = g(X,Y,a), (15) 
at 

where we have used (6 X ) = (6 y ) = (6") = and defined the second moment Q* = ((5 X ) 2 ). We 
also defined Q? = ((6 y ) 2 ), a 2 = <(#p 2 >, and the shape factors = (6 x 6 y ), = (6 x 6 a ) and 
W = (8 y j5"). The evolution equations for the second moments are found from the first-order 
expansion 6 X - = xj - X, so that Q v = 2(8 X 8 X ) and S^' = (S x 5 y + 5 x 8 y ), were the dot stands for time 
derivative. 



(16) 



es x = m + ./;,>; + l -f xx [{6 x f - a x ] , 

K = 8^ + g y 6 y + g a d^ (17) 

Cr = - e [f x a x + f^], (18) 

a y = ilg^ + gyW + gjr], d9) 

= i [/^ + j^] + + ^ + gX\ (20) 

t » = l[f x r a +f y I? a ], (21) 

27" = ^+ ft ,2>' fl +^cr 2 . (22) 

The system of Eq. (fT4l)-(fl~5l) together with Eq. (fT8l)-(l22l) forms a closed set of differential 
equations for the average collective variables X(t) and Y(t): 



eX 
Y 
eCl x 



-X 3 + (1 + b)X 2 -(b + 3Q, X )X + (1 + b)Q x + d - Y, 
In 



X - cY + a + A sin I — tj , 
2(-3X 2 + 2(1 +b)X-b- K)Q X - 2Z xy , 
2 [S^ - cQ? + I?' a ] , 

- [(-3X 2 + 2(1 + b)X - b - K)J?» - Q. y ] 



(23) 
(24) 
(25) 
(26) 
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+ Q x - clP' + S xfl , (27) 
e± xa = (-3X 2 + 2(\+b)X-b- K)L xa - Z yQ , (28) 

jya = jxa _ ^ya + ^2 (29) 

Numerical integration of this system allows us to obtain X(f), from which we can compute the 
spectral amplification factor 77. The value of rj obtained from the expansion is later compared with 
that obtained from the numerical integration of the Eqs. ©-© (see Figure [9] below). 

We can obtain another set of closed equations for X(t) and Y(i) if we perform an adiabatic 
elimination of the fluctuations, i.e., Cl x = Cl y = t, xy = t, xa = t, ya = 0, yielding to: 

2 

\xa _ °~ yya _ H(x)cr 2 ^xy _ H(x)cr 2 

CH(X)-V cff W-!' (cH(x)-l)^ 

qx _ °" ov _ H\x)a- 2 

(c//(x)-l) 2 ' ™-D 2 ' 

with H(x) = -3x 2 + 2(1 + b)x - b - K. Substituting Q x in Eqs. d23)-d2U), we find a closed form 
for the equations describing the evolution of the mean-field variables X(t) and Y(t): 



eX = -X 3 +(1 +b)X 2 



b + 



(cH(X) - l) 2 



X + (1+ ^ 2 1 + d-Y, (31) 
(cH(X) - l) 2 



Y = X-cY + a+Asm(^-^ (32) 

These equations provide a closed form for the nullclines of the global variables X and Y for the 
non-forcing case A = 0. They also reflect how diversity influences these variables. Figure [7] shows 
these nullclines Y\(X, cr) and Y 2 (X, a) of Eqs. (f3TT)-(f32l) respectively for a = 0.06 and different 
values of the diversity cr. It can be seen in the figure that the diversity changes the shape of the 
cubic nullcline Y\(X, cr) leading to a loss of stability of the fixed point of the system that, for a 
certain range of cr, becomes a limit cycle. To schematize the behavior of the global variables X 
and Y when the diversity changes, we show in Figure [8] the phase-portrait for different values of 
cr = 0.0, 0.5, 0.8, 1.0, 1.2 and 1.4 (corresponding to the values represented in Figure|7]). It can be 
seen that there is a range of cr for which the system exhibits a collective oscillatory behavior even 
in the absence of the weak external modulation. 

With the collective variable X(t) obtained from the adiabatic elimination we can estimate the 
spectral amplification factor rj. Figure [9] shows the results obtained from the numerical integration 
of Eqs. (|23l) -(l29l). together with the numerical simulation of the full system, Eqs. ©-© and the 
adiabatic approximation obtained from Eqs. (TJTI) and (1321) . It can be seen that our order parameter 
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FIG. 7: Nullclines of Eq. (f3TT)-(13^T> for different values of the diversity cr. Y\(X, <x) for <x: (a) 0.0 (red line), 
(b) 0.5 (green line), (c) 0.8 (blue line), (d) 1.0 (violet line), (e) 1.2 (cian line) and (f) 1.4 (yellow line). The 
nullcline Y%{X, a) of Eq. (1321) for a = 0.06 is represented with a black line. 

expansion is in good agreement with the numerical integration of the full system, even in the case 
in which the second moments are adiabatically eliminated. 



V. CONCLUSIONS 



We have studied the effect of the diversity in an ensemble of coupled neurons described by the 
FHN model. We have observed that an intermediate value of diversity can enhance the response of 
the system to an external periodic forcing. We have studied both electrical and chemical coupling 
schemes finding that the electrical coupling induces a larger response of the system to an external 
weak modulation, as well as existing for a wider range. In contrast, the chemical coupling scheme 
exhibits a smaller optimal amplitude and narrower range of response, however, for smaller values 
of the diversity. We have also found that the response of the system in the electrical coupling 
scheme strongly depends on the fraction of connected neurons in the system whereas it does not 
improve much above a small fraction of connected neurons in the chemical coupling scheme. 
We have also developed an order parameter expansion whose results are in good agreement with 
those obtained numerically for the electrically (diffusively) coupled FHN system. By an adiabatic 
elimination of the fluctuations we have found a closed form of the effective nullclines of the global 
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FIG. 8: Phase-space portrait for different values of cr. Grey line represents the evolution of {X, Y}. The 
black lines represent the nullcline Y2(X, a) for a - 0.06 and the cubic nullcline Yi(X, cr) for cr: (a) 0.0, (b) 
0.5, (c) 0.8, (d) 1.0, (e) 1.2 and (f) 1.4. 

collective variables of the system obtaining a simple expression of how the diversity influences the 
collective variables of the system. 

The microscopic mechanism leading the system to a resonant behavior with the external signal 
is as follows: in the homogeneous situation, where all the units are identical, the weak external 
modulation cannot induce any spike in the system. When the diversity increases, a fraction of the 
neurons enters into the oscillatory regime and, due to the interactions, pull the other neurons with 
them. This leads the system to an oscillatory collective behavior that follows the external signal. 
For larger values of the diversity, the fraction of neurons in the oscillatory regime decreases and 
the rest of neurons offer some resistance to being pulled by the oscillatory ones; thus, the system 
cannot respond to the external signal anymore. These results suggest that the diversity present 
in biological systems may have an important role in enhancing the response of the system to the 
detection of weak signals. 
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FIG. 9: Order parameter expansion versus numerical integration of the full system. An adiabatic approxi- 
mation is also included (see text). Two different periods of the external modulation have been considered: 
(a) T = 1.6 and (b) T - 1,11. Other parameter as in Figure [3] 
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